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ESTIMATION  OF  TEMPORALLY  AND  SPATIALLY 
VARYING  COEFFICIENTS  IN  MODELS  FOR 


INSECT  DISPERSAL 


H.  T.  Banks,  P.  L.  Daniel  Lamm  and  P.  M.  Kareiva 


ABSTRACT 

We  describe  techniques  for  estimating  temporally  and  spatially  dependent 
parameters  (including  coefficients  )  that  appear  in  general  transport 
models.  Convergence  properties  of  the  resulting  algorithms  are  given  and 
sample  computational  findings  with  test  examples  are  presented.  We  con¬ 
clude  with  a  summary  of  our  use  of  the  methods  analyzing  experiments  on  the 
movements  of  marked  flea  beetles  in  cultivated  arrays  of  the  cole  crop, 
collards  (Brass ica  oleraceae) . 


I.  INTRODUCTION 


Transport  equations  appropriately  model  numerous  biological  systems 
and  have  played  an  important  role  in  mathematical  biology  (cf  Rubinow 
[14],  Okubo  [12]).  Recently  such  transport  equations  have  been  increas¬ 
ingly  used  with  experimental  data,  ranging  from  developmental  biology 
(cf  Kauffman  et  al.  [10])  to  population  biology  (cf  Levin  [11]).  This 
connection  between  model  and  data  can  be  facilitated  by  parameter  iden¬ 
tification  schemes  (i.e.,  algorithms  for  generating  the  parameters  in 
a  model  that  provide  the  best  match  between  model  predictions  and  ob¬ 
served  data).  For  the  past  several  years  we  have  been  developing  and 
testing  spline-based  algorithms  for  identifying  parameters  in  transport 
equations  (cf  Banks,  Crowley  and  Kunisch  [1]).  Our  most  recent  efforts 
have  been  directed  towards  models  of  insect  dispersal  where, working 
with  insect  mark- recapture  data,  we  have  had  some  success  in  identifying 
constant  coefficients  and  spatially  varying  coefficients  in  transport 
equations  of  insect  movement  (Banks  and  Kareiva  [4]).  However,  our 
successful  efforts  involved  only  our  more  abbreviated  data  sets  (those 
spanning  one  or  two  days)  and  were  counterbalanced  by  failures  at  de¬ 
scribing  shifts  in  insect  distributions  over  the  course  of  three  days. 

We  hypothesized  that  our  models  failed  in  some  instances  because  they 
lacked  temporal  variation  in  the  parameters  reflecting  rates  of  insect 
movement  and  migration  (Banks  and  Kareiva  [4]). 

The  idea  of  temporally  varying  dispersal  rates  is  not  new  to  bio¬ 
logists.  In  general,  the  complication  of  temporal  variation  in  rate 
constants  or  model  parameters  is  ubiquitous  in  biology.  For  example, 
many  biological  processes  vary  with  time  because  the  environment  changes 


seasonally.  Biological  processes  may  also  vary  temporally  without  any 
driving  environmental  force;  a  familiar  example  is  the  inexorable 
senescence  of  organs  in  aging  Individuals.  Motivated  by  the  apparent 
temporal  variation  in  the  mobility  of  the  insects  in  our  experiments, 
and  the  recognition  that  temporal  variation  is  a  general  feature  of  bio¬ 
logical  processes,  we  have  made  efforts  to  extend  our  parameter  estimation 
algorithms  to  enable  us  to  treat  transport  equations  that  contain  time- 
varying  parameters. 

In  this  paper  we  describe  our  methods  for  identifying  both  spatially 
and  temporally  varying  coefficients  in  transport  equations.  Although 
the  particular  examples  with  which  we  illustrate  our  methods  here  repre¬ 
sent  scalar  equations,  our  techniques  are  readily  used  with  vector  systems. 
Beginning  with  a  basic  transport  model,  we  sketch  the  ideas  behind  our 
spline  based  techniques  and  summarize  the  main  convergence  results.  We 
then  report  on  tests  of  our  methods  on  numerical  examples  involving  co¬ 
efficients  that  simultaneously  vary  in  space  and  time.  Because  we  view 
our  methods  as  tools,  we  discuss  practical  computing  aspects  of  implementing 
these  methods.  We  then  summarize  results  obtained  when  we  applied  the 
methods  to  experimental  data  that  describes  the  dispersal  of  flea  beetles. 

It  was  this  data  that  initially  prompted  us  to  consider  models  with  time- 
varying  coefficients.  Since  this  paper  represents  the  culmination  of 
several  mathematical  and  biological  investigations,  we  have  sacrificed 
many  details  in  order  to  present  an  overview.  Readers  interested  mainly 
in  the  biological  issues  of  flea  beetle  dispersal  should  consult 
Kareiva  [7],  [8],  [9],  and  Banks  and  Kareiva  [4]  ;  for  more  details 
concerning  our  general  approximation  techniques,  including  convergence 
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proofs,  readers  may  consult  Banks  and  Daniel  [2]  and  Banks,  Daniel 
and  Kareiva  [3]  . 

II.  SPLINE  BASED  ESTIMATION  TECHNIQUES  FOR  TRANSPORT  EQUATIONS 

We  continue  here  our  investigations  (Banks  and  Kareiva  [4],  Sec.  2) 
of  general  transport  equations  which  describe  advection,  Fickian  diffusion, 
and  sink/source  mechanisms  containing  variable  and  possibly  unknown  rate 
parameters.  Specifically  we  consider  the  problem  of  estimating  the  un¬ 
known  functional  parameters  U,  V,  a,  0,  and  y  appearing  in  the  system 

(Vu)  =  (V  ~±)  +  olu  +  f,  t€  (0  ,T] ,  xe(0,l) 

(1)  u(t ,0)  =  u(t,l)  =  0 

u(0,*)  =  uQ(y)  . 

It  is  assumed  that  y  *  y(x)  and  that  V ,  V,  a,  and  0  are  functions  of 
(t,x)€(0,T] x(0,l) .  We  note  that  although  (1)  is  formulated  in  terms  of 
homogeneous  boundary  conditions,  our  ideas  are  sufficiently  general  to 
include  systems  with  nontrivial  boundary  conditions,  possibly  even  depen¬ 
dent  upon  unknown  parameters.  A  standard  transformation  (see  Banks  and 
Kareiva  [4])  may  be  used  to  reformulate  such  systems  as  special  cases  of 
(1)  and  we  shall  therefore  restrict  our  considerations  in  this  paper  to 
estimation  problems  for  (1). 

The  parameter  estimation  problem  associated  with  the  set  q=  (P,U,cx,0,y) 
of  unknown  parameters  appearing  in  (1)  is  the  following:  Given  spatially 

A 

distributed  observations  u^  £  1^(0, 1)  corresponding  to  times  t^,  i“l,...,L, 
find  q  €  Q  (for  Q  a  given  class  of  admissible  parameter  functions)  that 
minimizes  a  given  fit-to-data  criterion.  We  employ  here  a  least  squares 


criterion  (although  we  are  not  limited  only  to  functionals  of  this 
type)  of  the  form 

L  - 

(2)  J(q)  =  l  |u(t  ;q)  -  u. | 

i*l 

where  u  is  the  solution  to  (1)  corresponding  to  q€Q;  the  parameter  set  Q 

is  assumed  to  possess boundedness ,  compactness,  and  regularity  properties 

(these  assumptions  are  made  precise  in  Banks,  Daniel  and  Kareiva  [3])  . 

Due  to  the  infinite  dimensional  nature  of  both  the  state  u  and  the 

parameters  q  (and  the  computational  difficulties  associated  with  optima- 

zation  over  infinite  dimensional  spaces),  one  of  our  objectives  has  been 

to  develop  numerical  estimation  algorithms  based  on  approximation  schemes 

for  both  the  state  and  the  parameters.  We  describe  such  an  algorithm 

N 

here.  To  this  end  we  define  finite  dimensional  state  spaces  H  and 

M 

approximate  parameter  sets  Q  and  reformulate  the  parameter  estimation 
problem  in  a  computationally  more  tractable  finite-dimensional  setting. 

To  consider  state  approximations,  we  first  rewrite  (1)  in  weak  form 
in  the  state  space  1^(0, 1)  with  the  usual  inner  product  <*,*>.  For  a 
fixed  q«  (P,P,a,$ ,y)€Q,  this  weak  form  is  given  by 

<ut,<|>N  -  <l’u,D4>>  +  <PDu,D<t>  >  -  <au,<)>>  =  <f  ,4> >,  t£(0,T] 

(3) 

u(0)  *  Uq (y) 

which  solutions  u  of  (1)  must  satisfy  for  all  sufficiently  smooth  ♦  with 

<)> (0)  -  <J)  (1)  *  0.  Here  D  denotes  the  spatial  differentiation  operator. 

We  next  use  Galerkin  ideas  with  (3)  to  approximate  solutions  u  of  (1): 

N 

For  each  N,  we  construct  finite  dimensional  subspaces  H  of  l^fO.l), 

N  ,  N  N  1  N 

H  *  span  (B, , . . . ,  where  the  basis  elements  B,  are  linear  com- 
1  N+i  J 
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binations  of  standard  Cubic  B-spline  elements  (see  Schultz  [15],  for 

example)  defined  on  a  uniform  mesh  of  size  1/N.  The  linear  combinations 

N 

are  chosen  so  that  the  resulting  basis  elements  satisfy  the  homogeneous 
boundary  conditions  in  (1).  For  q  and  N  given,  we  then  seek  the  solution 


N  N  N 

»V>  ■  jh 


which  satisfies  the  Galerkin  equations 


<u*,*>  -  <VuN,DiJi>  +  <PDuN,Di|)>  -  <auN,i])>  =  <f  ,i»  t€(0,T], 

(4) 

<uN(0),4)>  =  <u0(y),i})>  , 

N  N 

for  all  40  .  Since  H  is  finite  dimensional,  it  is  easy  to  see  that 

(4)  is  an  ordinary  differential  equation  in  the  "Fourier  coefficients" 
N 

Wj(t);  we  defer  further  discussions  on  the  nature  of  the  approximating 

equations  to  a  later  section  where  we  provide  details  on  the  computer 

implementation  of  the  overall  algorithm. 

Our  approach  to  approximation  of  the  functional  parameters 

q  =  (P,l/,a,6,y) €  Q  is  similar  to  that  of  the  state  approximations  in 

M 

that  we  construct  finite  dimensional  sets  Q  that  "approximate"  (in  an 

appropriate  sense)  the  original  infinite  dimensional  parameter  set  Q. 

Our  theory  allows  for  a  wide  variety  of  classes  of  approximation 
M  M 

sets  Q  (we  need  not  have  Q  C  Q)  which  in  general  are  unrelated 
to  the  Ntfl  order  state  spaces  HN.  In  the  next  section  we  indicate 
how  we  have  successfully  used  cubic  and  linear  spline-based  approxi¬ 
mations  for  the  parameters;  a  more  precise  statement  of  these  general 
ideas  may  be  found  in  Banks  and  Daniel  [2]  and  Banks,  Daniel,  and 
Kareiva  [3]  . 

In  practice  we  combine  both  state  and  parameter  approximation 
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r 

ideas  to  formulate  an  approximate  estimation  problem:  For  given  levels 

-N  M 

of  approximation  N  and  M,  the  problem  consists  of  determining  q  £ Q 

n 

that  minimizes 

(5)  JN(q)  =  l  |uN(ti;q)  -  uj2 
i=l 

M  N 

over  q£Q  .  As  before,  u  is  the  solution  to  (4)  corresponding  to  q. 

Each  approximate  estimation  problem,  being  finite  dimensional,  enjoys 
significant  computational  advantages  over  our  original  estimation  prr  'em 
in  addition  we  have  the  following  convergence  results  (valid  under  j  i- 
cient  regularity  assumptions  on  f,  and  Q  —  see  Banks,  Daniel  and 
Kareiva  [3])  for  solutions  of  the  approximate  problems: 


-N 

Theorem:  For  N=l,2,...,  and  M=l,2,...,  there  exists  a  solution  q.. 

-  M 

N  M 

to  the  problem  of  minimizing  J  over  Q  .  In  addition,  there  exists 
a  convergent  subsequence  qu  -*•  q*  such  that  q*£Q  is  a  solution  to  the 

“j 

original  problem  of  minimizing  J  over  Q.  Finally,  we  obtain  state 

\  A  \ 

variable  convergence  u  (t;q_,  )  -*■  u(t;q*),  where  u  ,  u  are  solutions 

"j 

to  (4),  (3)  corresponding  to  qM  ,  q*  respectively, 

1 

The  proof  of  this  theorem,  which  relies  heavily  on  variational 
arguments  and  spline  estimates  (e.g.,  see  Schultz  [15]),  may  be  found 
for  a  general  transport  equation  in  a  two-dimensional  domain  in  Banks, 
Daniel  and  Kareiva  [3]  and  thus  will  not  be  presented  here.  Instead 
we  turn  to  a  discussion  of  computer  implementation  and  numerical  testing 
of  our  algorithms,  and  to  our  findings  upon  use  of  these  techniques  with 
flea  beetle  dispersal  data. 


III.  NUMERICAL  IMPLEMENTATION 

In  this  section  we  first  outline  some  essential  features  of  computer 
implementation  of  the  spline-based  approximation  ideas  introduced  above. 
Then,  before  examining  the  more  difficult  problems  associated  with  esti¬ 


mating  parameters  for  models  using  experimental  data,  we  shall  present 

our  findings  for  a  representative  test  example. 

For  given  values  of  M  and  N,  we  wish  to  consider  the  approximate 

N 

parameter  estimation  problem,  namely  that  of  minimizing  J  of  (5)  over 
M 

Q  .  In  developing  a  computational  package  to  effect  this  minimization, 

one  can  often  take  into  account  special  characteristics  of  the  state 
N  M 

space  H  and  the  parameter  set  Q  .  Turning  first  to  the  state  space 

approximations,  we  examine  equation  (4)  for  any  given  value  of  parameter 

q  =  (P,l/,ot,f3  ,y)  in  QM.  We  recall  that  uN(t)C  HN  may  be  written  uN(t)  = 

N+l  N  N  N 

.EjWj(t)Bj  (where  B_.  are  the  cubic  spline  basis  elements  defined  in  sec. 

N 

2);  substituting  this  expression  into  (4)  and  letting  v  =  B^,  i=l,...,N+l, 
we  obtain  the  system  of  ordinary  differential  equations  (ODE) in  w^(t) =col (w^(t)) 

j  J  s 

N  *N  N  N  N  N 

Q V(t)  =  (KN(t)  +  G  (t))w  (t)  +  F  (t) 

(6)  nN  Nfm  N 
Q  w  (0)  =  Wq . 

Here  QN,  KN(t)  and  GN(t)  are  (N+l)  *  (N-*-l)  matrices,  where  for  i ,  j=l , . . .  ,N+1 , 

nN  -  <bN  bN> 

Qi,j  "  j’  i  ’ 

=<rP(t,*)  DBN+ U(t,*)BN,DBN> 
i » j  3  3  1 


N  ^  N  N 

gi.j  ’  <aty‘i>  ■ 
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In  general,  these  matrices  are  sparse,  banded,  and  (in  the  case  of  Q^,  GN) 

N  N 

symmetric.  In  addition,  the  (N+l)-vectors  F  (t)  and  Wq  are  given  by 

FN(t)  =  col  «f(B,t,*),B^» 

N  N 

wQ  =  co1(<Uq,B^  > ) . 

We  remark  here  that,  with  the  exception  of  the  time-varying  nature  of  the 

N 

matrices  above,  the  ODE  (6)  in  w  is  essentially  the  same  as  the  corres¬ 
ponding  system  (equation  (12))  in  Banks  and  Kareiva  (4).  The  fundamental 
difference  between  the  two,  however,  involves  the  way  in  which  variable 
parameters  are  treated  in  each.  In  Banks  and  Kareiva  [4],  a  fixed  para- 
metri2ation  is  chosen  that  imposes  a  specific  a  priori  shape  for  the  co¬ 
efficients  (e.g.,  in  that  reference,  the  choice  of  parametrizations  for  U, 
namely  l/(x)  =  c(x~.  5) ,  constrains  l'  to  a  c  lass  of  af  f  ine  functions).  Here  we 

allow  the  iterative  estimation  scheme  itself  to  determine  the  shape  of 

M 

the  coefficients  by  searching  over  functions  in  a  class  Q  which  contains 

a  large  number  of  different  functional  shapes.  (For  the  special  case 

described  below,  we  use  linear  and  cubic  splines  to  generate  the  elements 
M  , 

of  Q  . )  We  shall  briefly  indicate  some  of  the  computational  aspects  of 
this  approach  and  refer  the  reader  to  Banks  and  Daniel  [2]  for  a  more 
thorough  treatment  of  these  ideas  (we  note  that  the  underlying  model 
equation  in  Banks  and  Daniel  [2]  differs  from  the  one  considered  here, 
but  that  the  approximation  ideas  for  variable  parameters  presented  there 
are  directly  applicable  here). 

To  explain  the  implementation  of  approximation  techniques  for  para¬ 
meters,  we  shall  simplify  our  presentation  by  focusing  on  the  parameter 
a  only  and  note  that  the  parameters  V,  l/,  0  and  y  are  treated  in  a  similar 
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manner.  To  illustrate  our  ideas,  we  assume,  in  the  examples  to  follow,  a 

separable  form  for  the  parameters  (i.e.,  a(t,x)  =  a^COo^Cx))  and  that 

'  M  M  M  M 

the  approximation  from  Q  for  a  takes  the  form  a  »  where 

(7)  ct“(t)  -  l  a”c”(t) 

k*l 

(8)  a2(x)  »  l  akCk(x). 

k=l 

M 

Here  Ck>  k*l,...,M,  are  linear  spline  basis  elements  based  on  a  mesh  size 

T/(M-1)  and  Ck«  k=l,...,M  are  standard  cubic  spline  (B-splines  —  see 

Schultz  [15])  elements  with  mesh  size  l/(M-3).  (We  note  that  the  orders 
of  the  approximation  for  and  need  not  be  related;  the  form  assumed 

here  is  only  for  ease  in  expostion.)  For  fixed  M,  N,  the  matrix  G^ft)  of  e- 
quation  (6)  has  entries 


where  ak,  ak,  k=l,...,M,  are  to  be  estimated.  We  note  that  the  inner 

products  <CkB^,B^>  may  be  computed  and  stored  prior  to  beginning  the 

M  ~M 

optimization  process,  so  that  as  the  ak>  ak  are  updated,  we  need  only 

N 

recombine  the  terms  in  (9)  to  obtain  the  current  value  of  G  (t) .  It  is 
easy  to  see  how  to  extend  these  ideas  to  the  matrices  involving  approxi¬ 
mations  for  V,  8  and  y. 

N  M 

To  minimize  J  over  Q  using  any  one  of  a  number  of  iterative  opti¬ 
mization  procedures,  it  is  necessary  to  select  start-up  values  for  the 

M  aM 

parameters  (i.e.,  for  the  coefficients  ak>  ak  in  expansions  such  as 
(7),  (8)).  For  the  examples  reported  here  we  used  the  minimization 
package  LMDIF1,  which  is  MINPACK’s  version  of  the  Levenberg  -  Marquardt 
modified  Newton  algorithm.  We  also  used  IMSL's  DGEAR  (stiff  ODE  solver) 


to  compute  the  solution  to  (6)  for  each  parameter  iterate,  and  the 

N 

IMSL  numerical  quadrature  package  DCADRE  to  compute  J  (the  norm 
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difference  between  u  and  sample  data  at  sample  times  t^).  All  the  com¬ 
putations  reported  on  here  were  executed  on  the  CDC  6600  at  Southern 
Methodist  University  and  the  IBM  370/158  at  Brown  University. 

We  now  consider  a  test  example  where  (synthetic)  data  generated  from 
a  known  solution  (see  Banks  and  Daniel  [2])  is  used  with  the  estimation 
algorithm;  by  comparing  the  parameters  identified  by  our  estimation  pro¬ 
cess  with  known  "true"  parameter  values,  we  are  able  to  evaluate  the  per¬ 
formance  of  our  algorithm.  In  this  test  example,  "true"  parameters  are 
given  by  P=20,  l/(t  ,x)  =  -100 (x- .  5)  (4-t ) ,  and  ct  =  6  =  Y  =  0,  while  the  "true" 
state  is  given  by  u(t  ,x)  *  -4Q0x(x-l)cos (t/2)  ,  0  £  t  £  2 ,  0  <_  x  <_  1 . 

Starting  from  initial  guesses  of  the  parameters,  we  apply  the  optimization 

11  N 

procedure  to  J  ,  where  the  distributed  sample  data  needed  to  compute  J  is 

A 

given  by  u^  =  u(t^ , • ) ,  t^  =  .25i,  i=l,...,8.  For  each  example  presented 
below  we  let  N=3.  (Further  details  on  the  implementation  of  test  examples 
may  be  found  in  Banks  and  Daniel  [2]  .) 

Example  1. 


(a)  We  estimate  V  only, 

holding  V  fixed 

at  its  true 

value.  We  allow  time 

varying  estimates  for  V 

,  letting  ^(t)  = 

k^ldkCk(t) 

where  M=4 

,  ,  M 

and  the  C. 

k 

are  defined  as  in  (7). 

A  time-varying  initial  guess 

of  £>=40- 

15t  yields 

the  following  results: 

d" 

do 

d* 

d* 

1 

2 

3 

4 

initial  guess 

40. 

30. 

20. 

10. 

estimated  values  (N 

=3,M=4)  20.0002 

20.0002 

20.0001 

20.0000 

"true"  values 

20. 

20. 

20. 

20. 

Execution  time  was  110  seconds  with  a  final  value  of 

JN(^)-6 

M 

.7  *  10~7. 

10 


V 


(b)  We  estimate  the  time-varying  part  of  V  only,  holding  other  parameters 

fixed  at  true  values;  i.e.,  we  let  V(t,x)  *  V1(t)V2(x)  an<*  estimate 

1  i  M 

l/^(t)  *  4-t  only.  As  before,  N=3,  M=4,  with  l^*(t)  =  ^E^v^C^(t).  Starting 
from  an  initial  guess  of  U  ~1  we  obtain  the  following  estimates: 


M 

M 

M 

M 

V1 

V2 

V3 

V4 

initial  guess 

1.0 

1.0 

1.0 

1.0 

estimated  values  (N=3,M=4) 

3.99994 

3.33329 

2.66664 

1.99999 

"true"  values 

4.0 

3.33333 

2.66667 

2.0 

The  execution  time  for  this  example  was  approximately  62.55  seconds  and 

N  N  — N  —6 

the  converged  value  of  J  was  given  by  J  (1/^)  =  1.2  x  10 

Example  2. 

We  repeat  the  above  examples  (same  start-up  values,  N=3,  M=4)  but 
consider  the  case  where  sample  data  is  available  only  at  discrete  spatial 
points  x^  in  (0,1).  The  least  squares  fit-to-data  criterion  in  this  case 
is  given  by 

(10)  J(q)  =  l  |u(t  x  ;q)  -  u  |2  , 
i,j  1  J 


with  the  corresponding  approximate  criterion 
(11)  JN(q)  =  l  |uN(t  ,x  ;q)  -  u  j2; 

i  *  j  J  2 

here  u,  u  are  solutions  to  (3),  (4)  respectively  and  £  R  are  observed 
data  points.  Repeating  example  1(a)  (where  we  estimate  P) ,  we  obtain  the 
following  converged  values. 


estimated  values  (N=3,M-4)  19.9998  19.9996  19.9997  19.9997 

~N  -3 

The  CP  time  here  was  475.33  sec.  while  J  (P^)  *  1.44  x  10  .  Repeating 
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example  1(b)  resulted  In  the  following  converged  values  of  with  a 

"“N  — n 

CPU  time  of  143.99  sec.  and  the  least  squares  criterion  value  J  (1/ .)  = 

M 

4.57  x  l(f 7. 


V1  2  3  4 

estimated  values  (N=3,  M=4)  3.99994  3.33329  2.66664  1.99999 

Example  3. 

Here  we  estimate  both  spatially  and  temporally  varying  parts  of 

l/(t,x)  *  { 4—  t ]  {-100(x-.5)  ]=  using  the  ’’pointwise"  criterion 

J  .  As  in  the  previous  examples  we  approximate  U  by  1/^,  and  also 

M  J- 

approximate  by  l^(x)  =  E^v^C^x)  where  is  defined  as  in  (8)  and 
M  is  fixed  at  M=4.  Since  is  'the  product  of  and  l/^,  at  least  one 
of  their  basis  coefficients  must  be  fixed  or  there  will  be  an  infinite 
number  of  possible  combinations  of  l/^1  and  that  still  yield  1^1^  =  ^. 
For  the  case  considered  here  we  fixed  v.^  and  v"  at  their  "true"  values 
throughout  the  entire  optimization  process.  Initial  guesses  for  the 
remaining  coefficients,  as  well  as  their  converged  values  (for  N=3) , 
are  given  below. 


~M 

~M 

M 

M 

M 

M 

V3 

V4 

V1 

V2 

V3 

V4 

initial  guess 

-5.0 

-15.0 

1.0 

1.0 

1.0 

1.0 

estimated  values  (N=3,M=4) 

-8.3333 

-24.9883 

4.0005 

3.3338 

2.6670 

2.0003 

"true"  values 

-8.3333 

-25.0 

4.0 

3.3333 

2.6667 

2.0 

Additional  numerical  findings  for  this  example,  as  well  as  numerous  other 
examples  (with  both  temporally  and  spatially  varying  parameters),  may  he  found 
in  Banks  and  Daniel  Q 2];  Section  41 

We  consider  in  the  next  section  the  estimation  of  parameters  using 
models  for  insect  dispersal  with  experimental  data  from  field  studies. 
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The  numerical  results  reported  below  were  generated  using  a  state  approxi¬ 


mation  index  of  N-32  and  parameter  approximation  index  of  M*4.  In  addition, 

~N 

for  convenience,  the  "pointwise"  least  squares  criterion  J  was  used  in  the 

N 

minimization  scheme  instead  of  J  since  sample  data  u^  was  readily  available 
in  discrete  form.  As  in  the  test  examples  above,  we  determined  that  the  use 


~N  N 

of  J  vs.  J  made  little  difference  in  the  outcome  when  one  is  estimating 
time-varying  parameters.  (We  repeated  several  of  the  computations  with  J 


N 


as  the  criterion,  where  the  distributed  data  u^  was  constructed  by  inter- 

A  A 

polating  linearly  between  data  points  u^,...,u  ;  here  p  is  the  number  of 

observed  spatial  locations.) 


IV.  TEMPORALLY  VARYING  INSECT  DISPERSAL 
Since  dispersal  is  important  to  the  population  dynamics  of  insects 
(cf  Stinner  £t  al.  [17],  Joyce  [6]),  mathematical  models  have  been  used 
extensively  to  investigate  the  consequences  of  differing  patterns  of  in¬ 
sect  movement  (Kareiva  [7],  (8],  [9];  Okubo  [12]).  The  development  cf 
accurate  models  may  be  especially  valuable  in  the  design  of  pest  manage¬ 
ment  schemes,  where  certain  cropping  arrangements  appear  to  promote  out¬ 
breaks  by  altering  pest  movements  (Risch  et  al.  [13]).  As  part  of  our 
long-term  study  of  insect  dispersal,  we  have  been  examining  population 
models  that  are  special  cases  of  equation  (1)  (cf  Kareiva  [7];  Banks  and 
Kareiva  [4]).  An  extensive  review  of  field  mark-recapture  data  has  indi¬ 
cated  that,  although  passive  diffusion  models  are  a  good  beginning,  models 
of  insect  dispersal  must  also  allow  for  rates  that  vary  in  space  and  time 
(Kareiva  [8]).  In  our  first  attempt  at  treating  spatially  varying  dispersal, 
we  analyzed  the  spread  of  marked  populations  of  flea  beetles  and  found  that 
transport  models  with  spatially  varying  convection  (or  oriented  movement) 
performed  significantly  better  than  passive  diffusion  models,  but  still 


""1 


not  as  well  as  we  would  like  (Banks  and  Kareiva  [4]).  We  hypothesized 
that  in  these  flea  beetle  mark-release  experiments,  coefficients  re¬ 
presenting  certain  mechanisms  were  time  dependent.  Here  we  use  our  new 
techniques  to  test  the  hypothesis  that  flc..  beetle  dispersal  and  migra¬ 
tion  varies  temporally  in  these  mark-recapture  experiments.  Beetle  move¬ 
ment  was  observed  in  1  m  x  80  m  cultivated  strips  of  col lard  patches. 

Since  there  were  no  food  plants  for  flea  beetles  outside  of  these  linear 
arrays,  their  local  movement  was  confined  to  short  hops  up  and  down  each 
linear  array.  This  justifies  the  following  model  as  a  reasonable  hypo¬ 
thesis: 

a  a2 

(12)  -  V(t)  -  <x(t)u,  t  £  (0,T],  x€(0,l), 

3x 

where  u  is  the  density  of  marked  beetles,  V  is  the  diffusion  coefficient, 
a  represents  population  "decay"  rate  due  to  death  and  emigration,  and  x 
is  the  rescaled  spatial  scale  (where  the  interval [0,1]  corresponds  to  an 
actual  interval  of  length  100m).  On  this  rescaled  interval  the  center  of 
each  experimental  array  corresponds  to  x=.5,  and  the  data  are  observations 
of  u  at  evenly  spaced  sample  stations  between  .2  and  .8.  Initial  data  is 
given  by  u(0,x)=0  for  all  x^.5  and  u(0,.5)  is  the  number  of  marked  beetles 
that  were  released  at  the  mid-point  of  each  array  to  begin  the  experiment. 

Since  less  than  1%  of  the  marked  beetles  reached  the  ends  of  the  cultivated 
arrays  (x *  .1  and  x *  .9),  and  since  there  was  no  food  to  support  beetles 
outside  the  arrays,  we  used  u(t,0)  *  u(t,l)  *  0  as  our  boundary  conditions. 
More  details  concerning  the  marking  and  recapture  procedure  and  flea  beetle 
biology  can  be  found  in  Kareiva  [7].  Our  subsequent  analyses  of  experimental 
data  and  equation  (12)  are  reported  in  terms  of  days  after  release.  In  fact, 
the  data  were  collected  in  terms  of  hours,  with  8  hours  (9  a.m.  -  5  p.m.) 
representing  a  beetle's  "activity  day."  Outside  of  this  activity  period 
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beetles  are  Inactive,  neither  feeding  nor  moving  (Kareiva,  personal  obser¬ 
vation).  To  explore  the  hypothesis  of  time-varying  dynamics,  we  systema¬ 
tically  examine  the  ability  of  equation  (12)  to  describe  experimental 
data,  using  different  combinations  of  V  constant,  V  varying  in  time,  a 
constant  and  a  varying  in  time. 


The  spline  algorithm  we  have  developed  estimates  the  coefficients  in 
equation  (12)  that  yield  the  "best  fit"  between  model  dynamics  and  experi¬ 
mental  data;  in  particular,  the  algorithm  minimizes  the  sums  of  the  squared 


differences  between  observed  beetle  di. sities  u^  and  predicted  beetle  den¬ 
sities  u  at  all  points  (tj,x^)  for  which  there  are  data.  As  a  standard 
against  which  to  compare  the  fit  of  the  model,  we  consider  the  null  hypo¬ 


thesis  that  beetle  density  is  simply  a  constant  normal  random  variable. 


Using  this  null  hypothesis,  we  have  evaluated  the  success  of  differing 
forms  of  equation  (12)  by  an  ad  hoc  modification  of  multiple  regression 
analyses  and  significance  tests  based  on  the  F-distribution.  The  basic 
idea  in  this  approach  is  the  "extra  sums  of  squares"  principle  (Draper 
and  Smith  [5]),  which  involves  examining  the  reduction  in  the  residual 
sum  of  squares  that  is  achieved  by  adding  parameters  to  the  model.  Al¬ 
though  for  such  non-linear  models  use  of  the  F-test  cannot  be  rigorously 
justified,  it  at  least  provides  a  guide  to  the  significance  of  our  findings. 

Following  the  protocol  of  traditional  regression  analysis,  we  define 
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analysis,  however,  our  model  does  not  hypothesize  that  beetle  density  is 
a  simple  function  of  selected  independent  variables.  Instead,  our  model 
is  a  partial  differential  equation,  the  solution  of  which  depends  on  ini¬ 
tial  data  and  boundary  conditions  (known  parameters)  and  on  (possibly 
time-varying)  coefficients.  These  time-varying  coefficients  are  the  un¬ 
known  parameters  we  are  seeking  to  identify.  In  our  statistical  analyses, 
we  calculated  the  degrees  of  freedom  (df)  in  each  model  identification 
run  as  simply  the  number  of  parameters  in  equation  (12)  that  were  esti¬ 
mated  by  our  algorithm.  For  example,  if  we  were  testing  a  model  with 
a  =  0  and  V  an  unknown  constant,  the  model  df  is  1.  If  we  sought  to 
estimate  a  constant  V  and  a  time-varying  a  that  was  to  be  represented 

M 

by  4  linear  splines,  then  our  model  df  is  1  (for  V)  plus  4  (for  the  a^, 

k»l,...,4  in  (7))  or  5  in  total.  Given  a  total  sums  of  squares  (here- 

~N 

after  TSSQ)  and  error  sums  of  squares  (J  of  Section  3;  hereafter 
referred  to  as  ERRSQ) ,  we  calculated  explained  sums  of  squares  (EXSSQ) 
as  TSSQ  -  ERRSQ.  In  traditional  regression  analyses  EXSSQ  must  always 
be  nonnegative.  When  testing  dynamic  models  such  as  equation  (12),  it 
is  possible,  however,  to  obtain  ERRSQ  >  TSSQ  and  consequently  EXSSQ  <  0. 
This  possibility  arises  because  a  dynamical  model  may  be  so  inappropriate 
that  treating  the  data  as  simple  normally  distributed  (white)  noise  provides  a 
better  description  than  "foolish"  dynamics.  Thus  the  examination  of 
particular  models  such  as  (12)  through  parameter  estimation  techniques 
is  not  an  exercise  in  simple  curve  fitting  with  guaranteed  success 
whenever  enough  parameters  are  included.  In  our  situation,  no  amount 
of  refining  of  parameter  fitting  can  resurrect  a  poorly  chosen  dynamic 


model. 
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Three  separate  dispersal  exoeriments  were  performed  with  the  flea 
beetle  Phyllotreta  striolata:  one  in  a  linear  array  with  3  m  between 
collard  patches,  one  with  6  m  between  collard  patches,  and  one  with  9  m 
between  patches.  At  all  three  interpatch  spacings,  our  EXSSQ  exceeded 
75%  of  the  TSSQ  when  both  P  and  a  were  allowed  to  vary  in  time  (Table  1). 
However,  a  good  fit  between  predicted  and  observed  data  could  in  fact, 
be  consistently  obtained  with  a  much  simpler  model  in  which  we  treated 
P  as  a  constant  but  allowed  a  to  vary  temporally.  In  this  case  we  always 
explained  a  significant  portion  of  the  TSSQ,  a  portion  which  was  only 
negligibly  less  than  the  portion  explained  when  P  also  varied  temporally 
(Table  2).  Without  temporal  variation  in  either  P  or  a,  solutions  of 
equation  (12)  could  provide  a  good  fit  with  data  for  only  one  day  (Figure  1) . 

An  interesting  methodological  point  concerning  the  search  for  time- 
varying  parameters  is  the  importance  of  initial  guesses.  In  Table  3  we 
see  that  using  the  same  experimental  data  and  the  same  model  but  differing 
initial  guesses,  we  obtained  dramatically  different  parameter  values  and 
percentages  of  TSSQ  explained.  Our  "good"  first  guesses  (those  eventually 
leading  to  large  EXSSQ)  were  the  constant  a's  or  P's  identified  by  our 
spline-based  algorithm  under  the  assumption  of  no  temporal  variation  in 
parameter  values.  These  estimates  of  constant  a  and  P  then  could  be 
effectively  used  as  starting  points  for  our  search  of  temporally  varying 
parameters. 

Not  only  was  a  consistently  identified  as  significantly  varying  in 
time,  the  quantitative  nature  of  that  temporal  variation  was  remarkably 
Identical  among  different  experiments  (Figure  2).  We  interpret  Figure  2 
to  represent  the  following  general  biological  phenomenon:  Early  after 
their  release  (l.e.,  in  the  first  day),  there  is  a  peak  in  beetle  emi¬ 
gration  (high  a)  due  to  the  disturbance  of  being  marked  and  handled; 

17 


after  one  day  of  natural  activity  in  the  field  this  disturbance  effect 
disappears  and  a  falls  to  its  "natural"  level.  Disturbance  following 
marking  and  release  is  so  common  in  insect  dispersal  experiments  that 
ecological  methodology  books  warn  against  it  as  an  experimental  artifact 
that  needs  to  be  controlled  for  or  taken  into  account  when  analyzing 
.data  (cf  .Southwood[16J  ;p.91)  .  Because  our  method  allows  us  to  explicitly 
identify  the  disturbance  effect,  in  studying  specific  models  we  have 
a  quantitative  technique  for  estimating  intrinsic  natural  diffusivity 
and  emigration  of  flea  beetles  after  disturbance  has  been  factored  out. 
We  are  now  applying  our  methods  to  insect  dispersal  systems  where  we 
suspect  that  day  to  day  variation  in  temperature  produces  marked  varia¬ 
tion  in  movement  rates.  Our  identification  algorithms  provide  a  tool 
for  describing  the  relationship  between  V  and  temperature  under  field 
circumstances  that  do  not  permit  the  experimental  control  or  manipu¬ 
lation  of  temperature. 
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TABLE  1.  Least  Squares  Analysis  of  the  Ability  of  Differing  Transport 
Models  to  Describe  Beetle  Dispersal  Data. 

DESCRIPTION  OF  DATA 

3  m  spacing:  9  points  in  space  sampled  1  and  3  days  after  beetles 
were  released,  TSSQ  *  76.3 

6  m  spacing:  9  points  in  space  sampled  1  and  3  days  after  beetles 
were  released,  TSSQ  ■  53.5 

9  m  spacing:  7  points  in  space  sampled  1  and  3  days  after  beetles 
wer^  released,  TSSQ  =  109.1 

_  PERFORMANCE  OF  VARIOUS  MODELS _ _ 


%  of  TSSQ 

Data 

explained  by 

F-scatistic  and 

Model 

Set 

model 

significance  level 

3u  n  3^u 

=  p - otu 

3t  .  2 

3  m 

32.7% 

F2  15  =  3.64,  p  >  .05 

3x 

frrith  constant 

6  m 

NONE 

(model  ERSSQ >  TSSQ) 

3>  and  a) 

9  m 

73.4% 

F2  u  =  15.97,  p<  .005 

3u  I,  32 u  ,  v 

S'’?'*’' 

3  m 

83.9% 

F3  12  =  12.48,  p<  .005 

(with  constant 

6  m 

70.4  % 

F5,12  =  5*71’  P<  -01 

Si,  variable  a) 

9  m 

97.6% 

F^  9  =  68.7,  p  <  .005 

|H  -  p(t)  _  a(t)u 

c  3x 

3  m 

92.5% 

Fgj9  =  13.8,  p<  -005 

(with  variable 

6  m 

76.7  % 

?8  9  =  3.71,  p<  .05 

Si,  variable  a) 

9  m 

97.6% 

F8,15‘1?-2'P<-°05 
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TABLE  2.  Choosing  the  "Best"  Model  By  Examining  How  the  Percent  of 

TSSQ  Explained  Can  Be  Increased  as  the  Models  Gain  Complexity 
(i.e..  Number  of  Parameters) 


Moving  from  constant  a  and  D  (2  parameters) to  a  constant  V  but 
variable  a  model  (5  parameters) _ _ 


3  m  spacing: 


6  m  spacing: 


9  m  spacing: 


%  TSSQ  explained  goes  from 
32.7  to  83.9 


%  TSSQ  explained  goes  from 
0  to  70.4 


%  TSSQ  explained  goes  from 
73.4  to  97.6 


F3,12  *  12’7'  P  <  -°05 

significant  improvement 

F3,12  =  9- 51*  p  <  '°05 
significant  improvement 

F3  g  =  18.82,  p  <  .005 
significant  improvement 


B.  Moving  from  constant  V  and  variable  a  model  (5  parameters)  to 
variable  V  and  a.  model  (8  parameters) _ 


3  m  spacing: 


%  TSSQ  explained  goes  from 
83.9  to  92.5 


F3  9  =  3.43,  p  >  .05 
improvement  is  not  significant 


6  m  spacing: 


7c  TSSQ  explained  goes  from 
70.4  to  76.7 


F3  -  8.11,  P  >  .5 
improvement  is  not  significant 


9  m  spacing: 


there  is  no  change  in  %  no  improvement 

TSSQ  explained 
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TABLE  3.  Dependence  of  Performance  of  Algorithm  on  Initial  Guess 


Data  from  6  m  spacing  experiments  (see  Table  1) 
TSSQ  =  53.5 


MODEL 


3u 

3t 


.2 

=  V  ~  -  a (t)u 
3x 


%  of  TSSQ  explained  by  model 


EX1 :  V  =  .03101,  INITIAL  a(t) =  .20 

SEARCH  FOR  VARIABLE  a  NONE  (.model  ERSSQ  >  TSSQ) 


EX2:  V  =  .00046,  INITIAL  a(t)  =  .2567 

SEARCH  FOR  VARIABLE  a 


70.4% 


MODEL 

ft  ■  fk 

-  a(t)u 

EX3:* 

INITIAL  P(t) : 
INITIAL  a (t)  = 

.0296 

0 

.0302 

.00008 

.0000 

SEARCH  FOR  VARIABLE  V 

AND  a 

12% 

EX4 :  * 

INITIAL  P(t) : 
INITIAL  a(t): 

.0296 

.0097 

.0302 

.0263 

.00008 

.4834 

.0000 

.0058 

SEARCH  FOR  VARIABLE  V 

AND  a 

NONE  (model  ERSSQ  >  TSSQ) 

EX5 :  * 

INITIAL  P(t)  = 
INITIAL  a(t) : 

.00046 

.2559 

.2486 

.0056 

.00014 

SEARCH  FOR  VARIABLE  V 

AND  a 

76.7% 

* FUNCTIONAL  INITIAL  GUESSES  EXPRESSED  IN  TERMS  OF  SPLINE  REPRESENTATION 
COEFFICIENTS  -  SEE  SEC.  3 


2) 


FIGURE  LEGENDS 


Figure  1 


Figure  2 


Fitting  models  to  data  on  the  dispersal  of  marked  flea  beetles. 
Observed  average  densities  of  recaptured  beetles  in  linear  array 
containing  collard  patches  at  9  m  spacing  are  represented  by 
solid  dots;  the  dashed  line  represents  best-fit  equation  (12) 
with  constant  D  and  constant  a;  the  solid  line  represents 
best-fit  equation  (12)  with  constant  D  and  variable  a. 

Temporal  variation  in  beetle  disappearance  rate  (a  of  equation  12). 


BEETLE  DENSITY  (nos./m) 


10 


DAY  I 


POPULATION  DECAY  CONSTANT,  <*(t) 


05  I  15  2  25  3 


TIME  AFTER  RELEASE  (DAYS) 


FIGURE  2 
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